Extinction events and species lifetimes in a simple ecological model 
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A model for large-scale evolution recently introduced by Amaral and Meyer is studied analytically 
and numerically. Species are located at different trophic levels and become extinct if their prey 
becomes extinct. It is proved that this model is self-organized critical in the thermodynamic limit, 
with an exponent 2 characterizing the size distribution of extinction events. The lifetime distribution 
of species, cutoffs due to finite-size effects, and other quantities are evaluated. The relevance of this 
model to biological evolution is critically assessed. 



In complex systems, not everyone is equal. In a human 
society, different individuals fulfil different roles. Simi- 
larly, companies occupy different niches in an economic 
system, and species occupy different niches in an ecosys- 
tem. Among the rich variety of possible structures in 
these systems, power laws take a prominent place and 
characterize the size distribution of cities ||l|, incomes 
1^, and ecological extinction events ||]. Several models 
that lead to such a scaling scaling behavior have been in- 
troduced in the literature (for some recent examples see 
iQ-^), however, the picture is far from complete. 

While reality is certainly best described by a compli- 
cated web of interactions at all levels, the simplest type of 
models that incorporate qualitative differences between 
individuals (or companies or species) are hierarchical or 
layered models, where individuals in a given layer affect 
individuals in the neighboring layer. Such a structure can 
e.g., be found in ecosystems, where species are at differ- 
ent levels in a food chain, or in economic systems, where 
different types of producers are located at different places 
in a production chain. Recently, Amaral and Meyer |^] 
introduced a model for large-scale evolution that contains 
several trophic levels. Using computer simulations, they 
found a power-law size distribution of extinction events. 
It is the intention of this paper to prove that this model is 
indeed critical in the thermodynamic limit, and to eval- 
uate some of its properties analytically and numerically. 

The model is defined as follows: Species can occupy 
niches in a model ecosystem with L levels in the food 
chain, and N niches in each level. Species from the first 
level Z = do not depend on other species for their food, 
while species on the higher levels I feed each on k or less 
species in the level I — 1. Changes in the system occur 
due to two processes: (a) Creation of new species with a 
rate fj, in each empty niche. If the new species arises in 
a level I > 0, k prey species are chosen at random from 
the layer below. A species never changes its prey after 
this initial choice, (b) Extinction: At rate p, species in 
the first level I = become extinct. Any species in layer 
I = 1 and subsequently in higher levels, for which all 
preys have become extinct, also become extinct imme- 



diately. This rule leads to avalanches of extinction that 
may extend through several layers and will be shown be- 
low to obey a power-law size distribution. These rules 
are slightly different from 0, however, the results can 
expected to be the same. 

The dynamics of the model are characterized by slow 
driving (speciation), and by rare and fast avalanches (or 
extinction events). In this respect, the model is similar to 
models for sandpiles js), forest-fires P], and earthquakes 
Ip^ that are self-organized critical g with a power-law 
size distribution of relaxation events, if certain conditions 
are satisfied. 

Let us first discuss the case k — 1, where each species 
in layers I > feeds on one prey species only. In this 
case, each species is connected to exactly one species in 
layer I = via a food chain. Since several species can 
feed on the same prey, the structure of the ecosystem 
looks like a set of trees, each consisting of all species that 
are connected to the same bottom species. If a bottom 
species has existed for a long time, the tree connected to 
it extends through many layers and consists of a large 
number of species. If a bottom species is young, only few 
species of the lowest layers are connected to it, and the 
corresponding "tree" is small. When a bottom species 
becomes extinct, the whole tree of species connected to 
it becomes also extinct. Since each bottom species be- 
comes extinct with the same rate p, the size distribution 
of extinction events is identical to the size distribution of 
trees. 

After some time, the system can be expected to be in 
a stationary state where the speciation and extinction 
rates balance each other within each layer, leading to a 
constant mean species densities. Of course, for a finite 
system size N, there may exist considerable fiuctuations 
around the mean density. Let pi denote the species den- 
sity in layer I. The equation of motion for pi is 



dpi/dt = /i(l - pi) -ppi , 



(1) 



leading to a stationary density pi = p/{p + p) in each 
layer. The lifetime distribution of species is an exponen- 
tial, 
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PT{T)=peM-pT), (2) 

and is the same for each species. Let Si^\t) denote the 
number of species in layer / that are connected to bottom 
species i. Since each newly created species picks its prey 
species at random among the existing species in layer 
Z — 1 , the growth of si*-* {t) is given by 

ds'i'^/dt = pi)/pi-M-i = p4-i . (3) 

as long as bottom species i does not become extinct. The 
second identity holds in the stationary state. The size of 
a "tree", S'^''> = E; obeys consequently dS^'^l/dt = 
p5'(*) , leading to 

= 5W(0)exp(pi), (4) 

where the initial tree size is S'(*)(0) = 1, since only the 
bottom species is present (the age of a tree is measured 
from the creation of the bottom species on). 

The size distribution of trees, P{S), is related to the 
age distribution of trees, P{t), via 

P{S) = P{t)dt/dS . (5) 

Now, P{t) obeys the equation dP{t)/dt = —pP{t), or 
P{t) oc exp{~pt) , which can be combined with Eqs. (|^) 
and (0) to give 

P{S) oc . (6) 

Since the size distribution of trees is identical to the size 
distribution of extinction events, the exponent character- 
izing the extinction events is r = 2, in excellent agree- 
ment with the numerical findings by Amaral and Meyer 

0. 

There are two types of finite-size effects that modify 
this power law: First, there are effects due to the finite 
size of N, when the number of layers L is large. Since 
the height of a tree is proportional to its age, the number 
of trees extending up to layer I decrases as exp(— /), and 
the mean number of species within layer I that belong to 
the same tree is proportional to exp(Z). Thus, all species 
in layers / > In belong to the same tree and are simul- 
tanously destroyed. Consequently trees cannot become 
much higher that In layers, and there is a cutoff to the 
power law P{S) oc 5""^ at S^^ax NlnN. If the number 
of layers L is of the order InA^, the size of the largest 
extinction event is of the order NL, which is the total 
system size. Indeed, events of this size are reported in 
0, where N = 1000 and i = 6 ~ In TV. (These authors 
probably have chosen fc > 1.) 

Since real ecological systems have only a few trophic 
level, it is important to also study the case L < \nN, 
where the properties of the system do not depend on N 
and finite-size effects are due to the finite number of lay- 
ers, L. Since the height of trees grows linearly in time. 



while their size S grows exponentially, the typical size of 
the largest trees is now given by Smax oc expL. On the 
other hand, each extinction event destroys on an aver- 
age the same number of species in each layer, since the 
same number of species are created in each layer. Since 
each extinction event destroys one species in the bottom 
layer, the mean number of species destroyed per extinc- 
tion event is S — L. Together with the relation 




" SP{S)dS 



and the assumption that P{S) obeys a power law, this 
leads again to P{S) oc 5'"^ for S < Smax- 
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FIG. 1. The size distribution of extinction events for k — 1, 
p = 0.05, fi = 0.02, and L = 4 (solid line), L = 5 (dashed 
line), and L = 6 (dotted line). The straight line is a power 
law with the exponent t — 2. 

Figure 1 shows the size distribution of extinction events 
for different values of L. The size S is scaled by exp(i), 
and curves for different L collapse nicely. The number of 
layers in this simulation is not large enough to show the 
power law exponent r over a large scaling regime, but it 
is chosen such that it is close to the number of layers in 
real ecosystems. 

The results obtained so far hold for a broad range of 
values of the parameters fi and p. Of course, since fi and 
p are defined as rates, one must make sure in simulations 
with discrete time steps that they are small enough so 
that no artificial effects occur due to the parallel updat- 
ing of many sites. Also, there is assumed to be a time 
scale separation between extinction events which are fast 
on evolutionary time scale, and the creation rate of new 
species, which is much slower. Only in this case one can 
neglect the interference between new speciations and the 
extinction avalanches. 

Let us next discuss the properties of the model in the 
case A: > 1. Each newly created species feeds on k species 
in the layer below, and becomes extinct only after all k 
prey species have become extinct. Older species feed on 
less prey species and are therefore more likely to become 
extinct during a given time interval than younger species. 
Also, the prey of older species is in general older than the 
prey of younger species, and vanishes faster than the prey 
of younger species. In contrast, for k ^ 1 each species 
becomes extinct with the same probability, irrespective 
of its age. Choosing fc > 1 introduces correlations in 
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age and extinction probability between species and their 
preys that make analytic treatment harder. 

Figure 2 shows the size distibution of extinction events 
for k = 3, L ~ 6, and N = 999. As in a similar plot in 
, this size distribution appears to be a power law with 
extinction events up to the system size. In order to de- 
cide whether the system is indeed in a critical state, or 
whether the apparent power law is just due to a lucky 
choice of parameters, we have to discuss the system in 
the thermodynamic limit N, L oo with L < hiN. 



species for large / is obtained from Eqs. (0) and (||) by 
dropping the /-dependence in Pt and Pr- After a few 
steps, one obtains the following differential equation for 
Pr-- 



dPrir) 



[1 -TF.,(T)]i-i/fepr(T) 



-kdr/f. 



(9) 
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Integration of both sides, together with the condition 
Pt{0) = 1/T gives an implicit solution for Pt- For very 
large t, the first factor in the denominator on the left- 
hand side is close to 1 and can be neglected, leading to 



Pt{t) - exp(-A:T/r) 



for large t. On the other hand, for very large t the life 
time of a species is only limited by the lifetime of the bot- 
tom species to which it is connected. The bottom species 
is destroyed with a rate p, leading to k/T = p, or 



FIG. 2. The size distribution of extinction events for fc = 3, 
L = 6, N = 999, p = 0.05, and p = 0.01. The straight line is 
a power law fit with an exponent r — 2.009. 

Let P^\t) denote the lifetime distribution of newly 
created species in layer I. The following calculation re- 
lates p!ji\T) to P^ji~^\T). The preys of a newly created 
species in layer I are chosen at random from all the species 
in layer I — 1. Since species are created with a constant 
rate in each layer, the distribution of the remaining life 
time T of the prey species of a newly created species is 
given by Pr'~^^ (t), which is related to P^~^\t) via 



j^dTp:^-^\T) 

!^dT!:'dTPt'\T) 

(r)/T'( 



dTPr 



ii-i) 



i-i] 



(7) 



where T^'~^) is the mean lifetime distribution of species 
in layer I. The probability that a species created at time 
T = lives for a time T is identical to the probability 
that the last of its prey species becomes extinct at time 
T, leading to 



p!^\t) = 



1 - 

d_ 
df 



P!^'~^HT)d7 



k-1 



kp(i-^){T) 



-I k 



dTP^^-^\T) 



(8) 



T = k/p. 



(10) 



The asymptotic life time distribution of newly created 



Thus, after a few transient layers, the mean lifetime of 
newly created species saturates at k/p, which is k times 
the lifetime of species in the bottom layer. Figure 3 shows 
the lifetime distribution of species in the higher layers for 
different values oik, as obtained by numerically iterating 
the recursion relations Eqs. ^ and (H). One can see that 
the distribution becomes more peaked with increasing k. 

The fact that the life time distribution of species does 
not change any more after a few transient layers means 
that each extinction event destroys on an average the 
same number of species on each of the higher levels. Con- 
sequently the rate of species production is also the same 
at each of these levels. Let us now define a "tree" with 
index m to be the set of all species that would go extinct 
if bottom species m went extinct. Only species with one 
prey belong to those trees, and there are species with 
one prey that do not belong to a tree. The growth rate 
of such a tree must be proportional to its size, just as 
in the case k — 1, since in each layer (above the transi- 
tional layers) species with one prey are generated (from 
species with two preys) at the same rate. We can now re- 
peat the derivation of P{S) and Smax from above (case 
k — 1), and we find again t — 2 and Smax exp(L). 
The species density in levels above the transient levels is 
p = ^i/{^i+plk). 
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FIG. 3. The lifetime distribution of species in the higher 
layers for k = 1 (solid), fe = 2 (dotted), k = 3 (dashed), and 
= 4 (long-dashed). The time axis is scaled by 1/k. 

One of several unrealistic features of this model is that 
the number of species does not decrease at higher levels. 
However, the biomass must decrease exponentially with 
the level number, since not 100 percent of the biomass 
at a given level are consumed by the species in the next 
level, and since not all of the consumed prey mass is 
turned into predator mass. One might counter this ar- 
gument by saying that the dependencies between species 
are not limited to predator-prey relationships, and that 
extinction avalanches will not only pass from preys to 
predators, but also to many other species. This leads, 
however, to interaction loops instead of nice hierarchies, 
and the hierarchical model must be viewed as some crude 
mean-field-like approximation to the more complicated 
reality. 

If rule (b) is modified such that a species goes extinct 
as soon as one of its k prey species goes extinct, the 
number of species decreases exponentially with the level 
number: Each species is now connected to fc' bottom 
species, so that the death rate of species increases with 
layer number as A:' /p. If the speciation rate is constant in 
each layer, the species density decreases as while the 
size distribution of extinction events is still a power law 
with T = 2 in the thermodynamic limit. However, this 
power law cannot be seen for the system sizes used in the 
simulations. If the speciation rate is chosen to be propor- 
tional to the density of species in a layer, the density in 
the stationary state decreases as fc"'*^'"^-*/^, and the size 
distribution of extinction events has a cutoff after a few 
layers, even in the thermodynamic limit. 

To summarize, the model discussed in this paper is 
self-organized critical with a power-law size distribution 
of extinction events in the thermodynamic limit. Finite 
systems with only a few layers show this power law only if 
k is larger than 1 or 2, and a modified version of the model 
either does not show a power law in systems with few lay- 
ers, or is not critical at all. Thus, power law extinction 



events are not a generic feature of food-chain models in 
general, but occur only in some versions of these models. 
Also, a more detailed model that includes adaptation 
of species to their prey and that evaluates the transfer of 
resources from one layer to the next, was shown to be 
not critical. A model that is completely different from 
the one discussed in this paper, but equally simple was 
introduced some time ago by Bak and Sneppen Q . This 
model does not incorporate any layered structure, but it 
includes the fitness of species, and it gives a power law 
with an exponent different from 2. 

While the study of simple models like the ones men- 
tioned here is a necessary stage in the attempt to under- 
stand complex phenomena like large-scale evolution, all 
of them are unrealistic in many respects, and it can be 
doubted that they are capable of grasping all important 
features of evolutionary dynamics. Certainly, far more 
research is needed to gain a deeper understanding of the 
processes that lead the the observed patterns in the fos- 
sile record. 
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